function [WeiDist] = Principal_A2_n(Par1)
n = 50;
load matrix_time1_eh1_comp0.out
load stock_time1_eh1_comp0_anyhs0.out
load stock_time1_eh1_comp0_anyhs1.out
load stock_time1_eh1_comp0.out
load trans_time1_eh1_comp0.out
load W_time1_eh1_comp0.out

% Health transitions and unhealthy proportion
phl_data=matrix_time1_eh1_comp0(1,1);
plh_data=matrix_time1_eh1_comp0(1,2);
prob_l=matrix_time1_eh1_comp0(1,3);

w_mom=W_time1_eh1_comp0(1,:); % Read wage moments
meanwi_1=w_mom(1); meanwf_1=w_mom(2);   meanwi_2=w_mom(3);  meanwf_2=w_mom(4);  sdwi_1=w_mom(5); sdwf_1=w_mom(6);  sdwi_2=w_mom(7);   sdwf_2=w_mom(8);
p10wi_1=w_mom(9);  p10wf_1=w_mom(10); p10wi_2=w_mom(11); p10wf_2=w_mom(12); p25wi_1=w_mom(13); p25wf_1=w_mom(14); p25wi_2=w_mom(15); p25wf_2=w_mom(16); 
p50wi_1=w_mom(17); p50wf_1=w_mom(18); p50wi_2=w_mom(19); p50wf_2=w_mom(20); p75wi_1=w_mom(21); p75wf_1=w_mom(22); p75wi_2=w_mom(23); p75wf_2=w_mom(24); 
p90wi_1=w_mom(25); p90wf_1=w_mom(26); p90wi_2=w_mom(27); p90wf_2=w_mom(28); minwi_1=w_mom(29); minwf_1=w_mom(30); minwi_2=w_mom(31); minwf_2=w_mom(32); 
maxwi_1=w_mom(33); maxwf_1=w_mom(34); maxwi_2=w_mom(35); maxwf_2=w_mom(36);

wf_1=linspace(minwf_1,maxwf_1,n)';        % Wage in formal sector
wi_1=linspace(minwi_1,maxwi_1,n)';        % Wage in informal sector
wf_2=linspace(minwf_2,maxwf_2,n)';      % Wage in formal sector
wi_2=linspace(minwi_2,maxwi_2,n)';      % Wage in informal sector

wage_moments_data=[meanwi_1; meanwf_1; meanwi_2; meanwf_2; sdwi_1; sdwf_1; sdwi_2; sdwf_2;...
    p10wi_1; p10wf_1; p10wi_2; p10wf_2; p25wi_1; p25wf_1; p25wi_2; p25wf_2;...
    p50wi_1; p50wf_1; p50wi_2; p50wf_2; p75wi_1; p75wf_1; p75wi_2; p75wf_2;...
    p90wi_1; p90wf_1; p90wi_2; p90wf_2];

Dij=trans_time1_eh1_comp0(1,:)'; % Read transition moments

stock_mom_h=stock_time1_eh1_comp0_anyhs0(1,1:9); % Read stocks [FF_h;FI_h;FN_h;IF_h;NF_h;II_h;IN_h;NI_h;NN_h]
stock_h=stock_mom_h';

stock_mom_agg=stock_time1_eh1_comp0_anyhs0(1,1:9).*(1-prob_l)+stock_time1_eh1_comp0_anyhs1(1,1:9).*prob_l;
stock_agg=stock_mom_agg';

% Bootstrap variances
phl_mom_boot=var(matrix_time1_eh1_comp0(:,1));
plh_mom_boot=var(matrix_time1_eh1_comp0(:,2));

D_mom_boot=var(trans_time1_eh1_comp0)';
trans_data_boot=[ D_mom_boot(3); D_mom_boot(5); D_mom_boot(1); D_mom_boot(2); ...
    D_mom_boot(4); D_mom_boot(6); D_mom_boot(11); D_mom_boot(13);...
     D_mom_boot(9); D_mom_boot(10); D_mom_boot(12); D_mom_boot(14); ...
     D_mom_boot(15); D_mom_boot(16); D_mom_boot(7); D_mom_boot(8)];


stock_mom_boot_h=stock_time1_eh1_comp0_anyhs0(:,1:9);
HHstates_data_boot_h=var(stock_mom_boot_h)';

stock_mom_boot_agg=stock_time1_eh1_comp0_anyhs0(:,1:9).*(1-prob_l)+stock_time1_eh1_comp0_anyhs1(:,1:9).*prob_l;
HHstates_data_boot_agg=var(stock_mom_boot_agg)';
HHstates_data_boot=[HHstates_data_boot_h;HHstates_data_boot_agg]; % stock variances

w_mom_var=var(W_time1_eh1_comp0)';
meanwi_1v=w_mom_var(1); meanwf_1v=w_mom_var(2);   meanwi_2v=w_mom_var(3);  meanwf_2v=w_mom_var(4);  sdwi_1v=w_mom_var(5); sdwf_1v=w_mom_var(6);  sdwi_2v=w_mom_var(7);   sdwf_2v=w_mom_var(8);
p10wi_1v=w_mom_var(9);  p10wf_1v=w_mom_var(10); p10wi_2v=w_mom_var(11); p10wf_2v=w_mom_var(12); p25wi_1v=w_mom_var(13); p25wf_1v=w_mom_var(14); p25wi_2v=w_mom_var(15); p25wf_2v=w_mom_var(16); 
p50wi_1v=w_mom_var(17); p50wf_1v=w_mom_var(18); p50wi_2v=w_mom_var(19); p50wf_2v=w_mom_var(20); p75wi_1v=w_mom_var(21); p75wf_1v=w_mom_var(22); p75wi_2v=w_mom_var(23); p75wf_2v=w_mom_var(24); 
p90wi_1v=w_mom_var(25); p90wf_1v=w_mom_var(26); p90wi_2v=w_mom_var(27); p90wf_2v=w_mom_var(28); minwi_1v=w_mom_var(29); minwf_1v=w_mom_var(30); minwi_2v=w_mom_var(31); minwf_2v=w_mom_var(32); 
maxwi_1v=w_mom_var(33); maxwf_1v=w_mom_var(34); maxwi_2v=w_mom_var(35); maxwf_2v=w_mom_var(36);
wage_data_boot = [meanwi_1v; meanwf_1v; meanwi_2v; meanwf_2v; sdwi_1v; sdwf_1v; sdwi_2v; sdwf_2v;...
    p10wi_1v; p10wf_1v; p10wi_2v; p10wf_2v; p25wi_1v; p25wf_1v; p25wi_2v; p25wf_2v;...
    p50wi_1v; p50wf_1v; p50wi_2v; p50wf_2v; p75wi_1v; p75wf_1v; p75wi_2v; p75wf_2v;...
    p90wi_1v; p90wf_1v; p90wi_2v; p90wf_2v];

boot_vector = [trans_data_boot; wage_data_boot; HHstates_data_boot; phl_mom_boot; plh_mom_boot];


% Transitions
Dnf_1	 =	Dij(1);
Dni_1	 =	Dij(2);
Dfn_1	 =	Dij(3);
Dff_1	 =	0;
Dfi_1	 =	Dij(4);
Din_1	 =	Dij(5);
Dii_1	 =	0;
Dif_1  =	Dij(6);
Dni_df_1 = Dij(7);    % Transitions for spouse 2 (spouse) from not working to informal when spouse 1 (head) loses job in formal sector
Dni_di_1 = Dij(8);    % Transitions for spouse 2 (spouse) from not working to informal when spouse 1 (head) loses job in informal sector

Dnf_2	 =	Dij(9);
Dni_2	 =	Dij(10);
Dfn_2	 =	Dij(11);
Dff_2	 =	0;
Dfi_2	 =	Dij(12);
Din_2	 =	Dij(13);
Dii_2	 =	0;
Dif_2  =	Dij(14);
Dni_df_2 = Dij(15);
Dni_di_2 = Dij(16);

 
%HHstates_data_reduced=[stock_h;stock_l]; % [FF_h;FI_h;FN_h;IF_h;NF_h;II_h;IN_h;NI_h;NN_h;FX_l;IX_l;NF_l;NI_l;NN_l]  joint-job state
HHstates_data=[stock_h;stock_agg]; % [FF_h;FI_h;FN_h;IF_h;NF_h;II_h;IN_h;NI_h;NN_h;FF_agg;FI_agg;FN_agg;IF_agg;NF_agg;II_agg;IN_agg;NI_agg;NN_agg]  joint-job state

Par1(21) = 10^Par1(21);
Par1(22)=1000*Par1(22);
Par1(23)=1000*Par1(23);
Par1(24)=1000*Par1(24);
Par1(25)=1000*Par1(25);
Par1(26)=1000*Par1(26);
Par1(27)=1000*Par1(27);

transitions_data=[Dfn_1;Din_1;Dnf_1;Dni_1;Dfi_1;Dif_1;Dfn_2;Din_2;...
     Dnf_2;Dni_2;Dfi_2;Dif_2;Dni_df_2;Dni_di_2;Dni_df_1;Dni_di_1];
 % solve dist of wage offers, value functions and simulate trajectory (MSM)
[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l,iterVF]=solve_value_function_A1(Par1,...
    Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
if (iterVF > 400)
WeiDist = 10000000000;
else
[HHstates,wage_moments,transitions,phl,plh]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,...
    wi_1,ff_2,wf_2,fi_2,wi_2,Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l);


MatCov = zeros(64);
for (i=1:64)
    MatCov(i,i) = 1/(boot_vector(i));
end
MomSim = zeros(64,1);
MomData = zeros(64,1);
MomData = [transitions_data; wage_moments_data; HHstates_data; phl_data; plh_data];
MomSim = [transitions; wage_moments; HHstates; phl; plh];
WeiDist = 20*(transpose(MomData - MomSim))*MatCov*(MomData-MomSim)
end
end
